Friction contribution to water-bond breakage kinetics 
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Based on the trajectories of the separation between water molecule pairs from MD simulations, 
we investigate the bond breakage dynamics in bulk water. From the spectrum of mean first-passage 
times, the Fokker-Planck equation allows us to derive the diffusivity profile along the separation 
coordinate and thus to unambiguously disentangle the effects of free-energy and local friction on 
the separation kinetics. For tightly coordinated water the friction is six times higher than in bulk, 
which can be interpreted in terms of a dominant reaction path that involves additional orthogonal 
coordinates. 
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I. INTRODUCTION 

The unique properties of liquid water are relevant for 
a broad range of processes in biology, chemistry, and 
physics, as well as for technological applications 
A prominent goal of recent research has been to re- 
late macroscopic properties (among those the notable 
anomalies and singularities) to the microscopic struc- 
ture and thus to the hydrogen (H) bonding pattern be- 
tween individual water molecules [2] . This goal has only 
partly been achieved. Indeed, even for the most ele- 
mentary kinetic process of breaking a single H-bond be- 
tween two water molecules that are embedded in the 
bulk liquid matrix, various viewpoints exist: In an early 
application of transition path sampling, it was found 
that in roughly half of the cases of an H-bond break- 
ing event a new bond forms right afterwards [s^l, sup- 
porting Stillinger's switching-of-allegiance description of 
the local water dynamics [2]. In later simulation works, 
the water reorientation during this H-bond switching was 
shown to occur quite abruptly j4i], in line with the pro- 
nounced rotational-translational motion coupling of indi- 
vidual water molecules § . The non-exponential H-bond 
relaxation was shown to be due to a coupling of bond 
making/breaking dynamics and the relative diffusion of 
water pairs [61], but not related to the local environment 
of H-bond forming water molecules 01 , which is surpris- 
ing in light of the above mentioned H-bond switching 
scenario. Clearly, the H-bond dynamics is intimately re- 
lated to the kinetics of e.g. protein folding [sj or solute 
dissociation ^ , so clarifying the kinetics of the binding 
and unbinding of water molecules is without doubt of 
fundamental importance. 

The concept of diffusion along a reaction coordinate 
(RC) has been fruitful for gaining insight into the under- 
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FIG. 1. [Color online] (a) Simulation snapshot visualized us- 
ing VMD [TBI ]: the coordinate R in the enlarged section is 
defined as the radial separation between the oxygen atoms, 
(b) Typical time series of R, the magnification reveals fluc- 
tuations on the subpicosecond scale. A simulation movie is 
provided as supplementary material [l6| . (c) Illustrative typ- 
ical reaction path involving in addition to the separation R 
an orthogonal component R± (cf. text in Sec. IV C)) . 



lying mechanisms of high-dimensional dynamics as in the 
case of protein foldin g, f or which various approaches to 
identify suitable RCs Il0|. 11| and to locate or character- 
have been developed. Here, 
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ize transition states [l^ili 
we choose the separation between two water molecules 
as the naive RC and show that a consistent description 
of the dynamics along the separation coordinate can be 
obtained. In fact, our stochastic analysis in terms of the 
Fokker-Planck (FP) equation with coordinate-dependent 
free-energy and diffusivity allows us to quantify to which 
extent degrees of freedom that are orthogonal to our cho- 
sen RC are involved in the reaction. As a main result, we 
find the relative translational friction in the first coordi- 
nation shell to be more than six-fold increased compared 
to bulk water. Application of transition rate theory with- 
out taking this local friction change into account under- 
estimates typical bond breakage times by a factor of two. 
Our analysis is based on 10 ns long trajectories of the 
separation R between the oxygen atoms of water pairs 
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provided by molecular dynamics (MD) simulations of the 
standard three point charge water model SPC/E see 
Fig. [IK) and b) for a snapshot and an example trajectory; 
a simulation movie is provided as supplementary mate- 
rial [l^ . We do not check for the presence of H-bonds — 
which would introduce an element of arbitrariness due to 
the H-bond definition fl8| — but rather base our discus- 
sion solely on the separation R; strictly speaking we do 
not consider H-bond kinetics but more generally water- 
bond kinetics. 

The paper is organized as follows: We start by review- 
ing the FP equation in Sec. [ill on which our analysis is 
based. The simulation setup is described in Sec. IIII A[ 
details regarding the trajectory analysis are given in 
Sec. IIII BI Results are presented in Sec. lIVI and discussed 
in Sec. [V] A summary of our main findings is given in 
Sec. I VII while technical aspects are covered in the appen- 
dices. 



II. FOKKER-PLANCK EQUATION FOR 
RADIAL DYNAMICS 



In the overdamped limit, the FP equation in three di- 
mensions describes the time evolution of the probability 
density ^ of observing a vectorial separation r at time t, 



d_ 

dt 



-^{r.t) = -V • J(r), 



(1) 



where the probability flux density 
J(r) = -*(r,t)V3D(r) • VU{r) - ^3D(r) • V*(r,i), 

has two contributions: (i) the overdamped motion due to 
an (effective) potential U and (ii) diffusion with a (pos- 
sibly) position dependent diffusivity tensor D 3D . Using 



the Einstein relation D ^-q = k^T fi connecting mo- 
bility and diffusivity, Eq. [T] can be rewritten as 



*(r,i) = V • (e-^^M^3D(r-) • V (e'^^M*(r,t) 



d_ 

dt 

(2) 

where (3 = l/{k^T) denotes the inverse thermal energy. 
Within this paper, we concentrate on the relative dynam- 
ics of two water molecules along their radial distance R. 
Since the diffusion tensor 



£>3D = I?0 



(3) 



D 



remains diagonal when introducing spherical coordinates 
(i?, 8, $) and the effective inter-molecular potential U 
due to symmetry depends on R only, the angular coordi- 
nates can be integrated out [l^. The time evolution of 
the radial probability distribution 



f2,TT fTT 

P{R,t)= d$ / de sinei?2^'(i?,e,$, 
Jo Jo 



t), (4) 



specifying the probability of finding a radial distance R 
at time t, is described by the simpler equation 



d_ 

dt 



PjR) 

i?2 



(5) 

where the pair radial diffusivity D may depend on i?. It 
is useful to absorb the factors R^ in Eq. [5] by defining a 
free-energy F = U — 2kBTlogR [l^ to recover the usual 
form of the one-dimensional FP equation [2^ [2l| 



(6) 



The free-energy F{R) = -fceT log (P(i?)) is obtained 
by Boltzmann inversion of the equilibrium probability 
(P(i?)). Determining D{R) is more subtle: Different 
procedures have been proposed in the context of pro- 
tein folding [23 - |25j or interfacial water diffusion f26l-[28|. 
Here, we obtain D{R) directly from measured mean first- 
passage times (MFPTs); for diffusive dynamics described 
by Eq. [6l the MFPT rfp of first reaching a separation i?t 
when starting off from R is given by [29| 



fit 



n,iR,R.) = J^ dP'^^^, 



^PF(R') /.fi' 



- di?"e-^^(«"), (7) 

I "'fimin 



assuming a reflective (zero-fiux) boundary condition at 
-Rmin < R < Rt- By differentiation, one readily gets 



D{R) 



dTipiR,Rt)/dR 



di?'e-^^(«') 



(8) 



Extracting MFPT curves rfp from simulation data thus 
allows to determine the separation dependent diffusivity 
D{R) governing the dynamics in the free-energy land- 
scape F{R); resulting diffusivity profiles are presented in 
SecHVl 



III. METHODS 
A. Simulation Setup 

MD simulations of the SPC/E water model are 
performed with the Gromacs simulation package {30| . 
Systems consisting of 895 and 2180 water molecules are 
simulated in a cubic box with periodic boundary condi- 
tions. At T = 300 K this corresponds to box sizes of 
roughly 3.0 x 3.0 x 3.0 nm^ and 4.0 x 4.0 x 4.0 nm^. Most 
simulations are performed using 895 molecules; results 
for target separations Rt = 1.9 nm stem from simula- 
tions involving 2180 molecules. We perform simulations 
at temperatures of T = 280, 300, 320 and 340 K for the 
small system and at T = 300 K for the large system at a 
pressure of P = 1 bar. At each temperature the system 
is equilibrated first in an NVT ensemble (constant par- 
ticle number, volume and temperature) for 100 ps and 
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then in an NPT ensemble (constant particle number, 
pressure and temperature) for 100 ps. Production runs 
are performed subsequently for t = 10 ns and configura- 
tions are saved each 10 fs for the small system and each 
100 fs for the large system. A Berendsen weak coupling 
thermostat and barostat jsij with a relaxation time of 
T = 1.0 ps is used for temperature and pressure con- 
trol. All non bonded interactions are cut off at a radius 
of Rc = 0.9 nm. Long-range electrostatic interactions 
are treated by the particle mesh Ewald method [s^ with 
tinfoil boundary conditions. An analytic long-range cor- 
rection for the Lennard- Jones interaction is applied to 
energy and pressure. 

The simulation movie available as supplementary ma- 
terial [iB] visualizes the relative dynamics of a chosen 
water pair in the MD simulation (T = 300 K) and was 
created using VMD [ll|. 



B. Molecular Dynamics Data Analysis 

The diffusion constant -DH2O of a single water molecule 
is determined from the long-time limit of the single- 
molecule mean squared displacement (MSD), i^HaO = 
limt_,oo {{r{t) - r(0))2) /(6t), as detailed in App. Ill 

The relative dynamics of all pairs of water molecules 
within the 10 ns long MD trajectory are resolved using 
their minimal image distance with a spatial resolution of 
Ai? = 0.002 nm and a temporal resolution of 5t — 20 fs. 
All paths starting within a distance /S.R/2 from R and 
crossing Rt — AR/2 at a time ifp later contribute to the 
mean first-passage time (MFPT) rfp = (<tp) and to the 
first-passage time (FPT) distribution /tp. Due to the 
periodicity of the system, the relative dynamics along the 
coordinate R is only meaningful for R < L/2 with box 
size L; we therefore only consider target distances Rt < 
L/2. Note that the absolute values of the MFPT curves 
T{p{R, Rt) sensibly depend on the time resolution St of 
the underlying trajectory; this sensitivity is discussed in 
App. E 

The derivative dT{p{R, Rt)/dR in Eq. [8]is determined 
by fitting a straight line to Tfp within a region of width 
0.032 nm around R (corresponding to 17 data points). 
The width of the region was empirically found to smooth 
out the statistical noise in the MFPT curves without hid- 
ing the relevant variations of the diffusivity. The integral 
in Eq. [8] is evaluated numerically, the equilibrium distri- 
bution {P{R)) is linearly interpolated and the reflective 
boundary set to iimin = 0.235 nm. Applying the same 
kind of procedure based on simulations of 2180 molecules 
in a cubic box of edge length L ss 4 nm allows to consider 
targets Rt up to 1.9 nm without introducing artifacts due 
to the periodicity of the simulation box and thus resolv- 
ing the diffusivity D over a larger range of separations 
R; finite size effects in the diffusivity profiles were not 
observed. 

We thoroughly checked that for all numerical steps of 
the data analysis, varying the spatial and temporal reso- 
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FIG. 2. [Color online] (a) Pair correlation function goo from 
MD simulations for various temperatures, (b) MFPT curves 
Tfp from MD data for T = 300 K and several target separa- 
tions Rt. (c) Diffusivity profiles D at 300 K from the dis- 
tribution in (a), the MFPTs in (b) and Eqs. [S] and [5] (same 
colors as in (b)). (d) Diffusivity profiles rescaled by the bulk 
difi^usivity 2 -DhsO for various temperatures. 



lutions as well as the position of the reflective boundary 
i?min had no significant impact on the resulting diffusiv- 
ity profiles. 



IV. RESULTS 



Mean First-Passage Times and Diffusivity 
Profiles 



Fig. (2^) shows the pair-correlation function goo{R) for 
different temperatures, the maxima indicating the po- 
sitions of the respective coordination shells. The free- 
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energy 



F{R) = -2kBT log R-ksT log {goo{R)), (9) 



exhibits a barrier of about Ik^T for crossing from the 
first to the second coordination shell as seen in Fig. [5^). 
MFPT curves rfp extracted from the simulation data for 
targets Rt ranging from 0.4 to 1.4 nm for T = 300 K are 
shown in Fig. [5]d) . They are converted, using Eq. |S1 into 
diffusivity profiles D{R) shown in Fig. [2]:); details of the 
trajectory analysis are given in Sec. lIIIBI There is rather 
good agreement between the curves for different target 
separations Rt , which is strictly expected only for a pure 
Markovian process as described by a one-dimensional FP 
equation. As will be discussed later on, this suggests that 
water bond breakage, defined as the passage from the first 
to the second coordination shell, is to a good approxi- 
mation Markovian. The deviations seen when R ^ Rt 
are expected, since on the short spatial scales associated 
with those first-passage events water motion is not dif- 
fusive; in fact, the crossover between ballistic and dif- 
fusive motion of single water molecules occurs at length 
scales of around 0.1 nm (cf. Fig. lAll in App. For 
increasing separation all curves saturate at a value equal 
to twice the diffusion constant of a single water molecule, 
D{R ^ oo) = 2Z)h2 ~ 5.1 nm-^/ns (denoted by a bro- 
ken line), as expected. 

As our main finding, the diffusivity profile exhibits a 
pronounced drop within the first coordination shell and 
reaches a minimum value of _D « 0.79 nm^/ns about 
six times smaller than in bulk, while factors ^ 2 were 
previously observed in simpler systems [33| . The thin, 
black line in Fig. [5J:) was obtained by evaluating MFPTs 
to a target separation Rt = 1.9 nm based on simulation 
data of the larger box with edge length L « 4.0 nm. 

Diffusivity profiles corresponding to distinct target 
separations i?t deviate from each other in two respects: 
(i) non-Markovian dynamics on short time and length 
scales lead to modifications for |it! — J?t| ^ 0.25 nm 
as discussed above, and (ii) the statistical uncertainty 
increases with increasing |_R — i?t| due to a decreas- 
ing number of recorded transition events contributing to 
the corresponding MFPTs. Smooth and reliable diffu- 
sivity profiles are thus obtained by joining the regions 
i?t— 0.45 nm < R < i?t— 0.35 nm of the diffusivity profiles 
corresponding to targets Rt = 0.6, 0.7, . . . , 1.4 nm. The 
resulting diffusivity profiles rescaled by twice the bulk dif- 
fusion constant Dn^o are shown in Fig. [5Jl) for various 
temperatures; the values of D^i^o at different temper- 
atures are in good agreement with previous simulation 
estimates and experiments as seen in Table HI 

Interestingly, deviations over a temperature span of 
80 K are very small; the main features of the profile, in- 
cluding the six-fold decrease within the first coordination 



TABLE I. Temperature dependence of the diffusion coefficient 
-DH2O of a single water molecule in bulk water. Simulation 
results for the SPC/E water model obtained by evaluation 
of the long time MSD (cf. App. [XJ are compared to results 
from previous simulation studies and to experimental findings 
(both with references). 



T [K] 


Db_20 
Simulations (SPC/E) 


[nm'^/ns] 

Experiments 


278 




1.313 [34] 


280 


1.60 ± 0.02 


1.44 [35] 


298 


2.75 [36], 2.70 [37] 


2.22-2.61 [34> 35, 38] 


300 


2.55 ± 0.05 




318 




3.575 [34] 


320 


3.70 ± 0.05 




340 


5.08 ± 0.05 




360 


6.60 ± 0.05 




shell. 


are accurately described by the heuristic formula 



D{R) « 2Dn,o [10-76 - 0.68 e'^^/^ _ 1 (27-50i?)^ 
+ 10 tanh (50(1 - AR)) - 0.34 tanh (l3.2 - 40^) 



-HO.ltanh [4.1 - lORj j, 

(10) 

where R = R/nm; Eq. [10] is shown as thin black line 
in Fig. [2JI). From the Arrhenius-like temperature de- 
pendence of the bulk diffusion coefficient (cf. Fig. IA2I in 
App. 1X1) it follows that the entire diffusivity profile obeys 
an Arrhenius law. 



B. Maxima in the MFPT-Curves at Small 
Separations 

According to Eq. [7] the MFPT-curve Tfp{R,Rt) is a 
strictly decreasing function of R; in contrast, as can be 
seen in the left panel of Fig. [3l which shows a close- 
up of the MFPTs of Fig. at small separations, the 
MFPT curves obtained from MD simulation data show a 
maximum at separations R « 0.26 nm. Since according 
to Eq. [8] a vanishing / positive slope of a MFPT-curve 
implies a diverging / negative diffusivity, the concept of 
Markovian dynamics obviously breaks down at such small 
separations. The diffusivity profiles in Fig. [21;) and d) are 
therefore only resolved for separations R > 0.265 nm. 

Though being counterintuitive at first sight, these 
maxima in the MFPTs can easily be understood by con- 
sidering the average oxygen-oxygen separation {R{t))j^^ 
of an ensemble of water pairs starting with defined ini- 
tial separation i?o at time t ~ 0. SPC/E- water molecules 
interact via Coulomb and via Lennard- Jones (LJ) inter- 
actions: for small separations the repulsive part of the 
LJ-potential significantly contributes to the total energy 
of a water pair (for R — 0.25 nm and T = 300 K: 
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FIG. 3. [Color online] Left: Enlarged view of the MFPT curves rfp from MD simulation data at small separations (same data 
as in Fig. Hja)). Right: Average oxygen-oxygen separations {R{t))j^ for water pairs with defined initial separation Ro at time 
t = 0. 



Uhj ~ 13.5 k^T). The corresponding water pair is thus 
expected to be quickly driven apart due to the repulsive 
LJ-force. The right panel of Fig. [3] indeed reveals that 
the average distance between water molecules starting 
at separations i?o ^ 0.25 nm increases strongly within 
fractions of picoseconds. The oscillations seen in the 
trace for Rq « 0.245 nm nicely match the time scale 
of inter-oxygen vibrations, which was found to be on 
the order of O.I — 0.2 ps Due to the repulsive LJ- 
interaction for R < 0.25 nm, the average separations 
of water pairs starting out at i?o ~ 0.245 nm and at 
i?o ~ 0.282 nm are very similar on time scales i ^ I ps 
and exceed the average separation of pairs starting in the 
range 0.255 nm < Rq < 0.275 nm; the maxima observed 
in the MFPT curves of Fig. [3] are a direct consequence 
of this mutual water repulsion and the induced under- 
damped motion. 



The accuracy of the FP approach involving the dif- 
fusivity profile is demonstrated when comparing first- 
passage time (FPT) distributions: Fig. ^p) contrasts the 
FPT histogram from MD data for Rt = 0.4 nm with 
FPT distributions from the numerical solution of the FP 
equation (numerical details are given in App. [C]), again 
using the flat diffusivity 2Dh2 and the actual diffusivity 
profile D{R). Only the FP approach including the D{R) 
profile correctly reproduces the entire FPT distribution 
from MD simulations and in particular also the exponen- 
tial tail of the distribution, as shown by the plot using 
the logarithmic scale on the right. Systematic discrepan- 
cies are observed on short time scales < 1 ps where the 
MD data show more "fast" transitions than the FP de- 
scription. These effects are caused by ballistic motion of 
water molecules and cannot be captured by a Markovian 
description. FPT distributions corresponding to other 
target separations are shown in Fig. ICIl in App. [C] 



V. DISCUSSION 

A. Fokker-Planck Kinetics with and without 
Diffusivity Profile 

To what extent is this local friction increase relevant 
for the water-bond breakage kinetics? To quantify the 
relevance of the change in local friction, we compare in 
Fig. 2^) MFPT curves from MD data already shown in 
Fig-llb) {colored lines) with analytical predictions result- 
ing from Eq. [7] using the diffusivity profiles D{R) {solid 
lines) shown in Fig. (^t) as well as calculations employ- 
ing a constant diffusivity D — 2_Dh2 {broken lines). The 
solid lines by construction match the MD data nicely, the 
vertical shift being caused by the 20 fs time discretization 
of the underlying MD trajectory, while the analytical pre- 
dictions are calculated in the continuum (cf. App. iBt. It 
is seen that the assumption of a constant diffusivity leads 
to a considerable underestimate of the MFPTs. The time 
to reach the target separation Rt = 0.4 nm from the first 
coordination shell {R < 0.28 nm) is underestimated by a 
factor of roughly one half. 



B. Testing the Quality of the Reaction Coordinate 

Although our procedure does not strictly depend on 
the fact that the separation i? is a "good" RC, the whole 
mapping on a one-dimensional FP equation is certainly 
more meaningful if it is. To check the quality of our 
RC, we divide R into a bound region A for R < Ra — 
0.275 nm, an unbound region B for R > Rb — 0.47 nm 
and the intermediate region for 0.275 nm < R < 0.47 nm 
which roughly encompasses the free-energy barrier (see 
Fig. [S^)). For a diffusive process described by the FP 
equation (Eq. [5]), the committor Trx{R) specifying the 
probability of first reaching region X G {A, B} when 
starting from i? is a solution of the stationary backward 
FP equation [U 



dR 



D{R) 



OMR) 
dR 



0. 



(11) 



The committor in addition fulfills boundary conditions 
7rx(i?x) = 1 and 7rx(i?Y) = 0, with Y^X. The solutions 
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FIG. 4. [Color online] (a) MFPTs from MD simulations {solid 
color lines, same data as in Fig. ^jp), and from the FP de- 
scription of Eq. [7] with constant diffusivity 2Dh2 (dashed) 
and with diffusivity profile D[R) from Fig. [2J;) [solid black 
lines) for several target separations Rt- (b) FPT distribution 
/fp to reach a separation Rt = 0.4 nm for water pairs starting 
within the first coordination shell at -R = 0.281 ± 0.001 nm. 
Histograms from MD simulations at T = 300 K (blue dots) 
are compared to the numerical solution of Eq.[5]with constant 
diffusivity 2_Dh20 (dashed black line) and with diffusivity pro- 
file D(R) from Fig. [2};) (solid black line). Data are shown on 
both linear and logarithmic vertical scales, vertical arrows 
indicate the mean of the corresponding distributions. FPT 
distributions for other target separations -Rt and numerical 
details are found in App. [Cl 



are given by 

I pR qPF{r') 
with the normalization factor 

pRb q0F(R') 



(12) 



(13) 



In Fig. \5jp) we show the committors 7ta{R) and tib{R), 
which specify the probability of first reaching region A 
and B, respectively, from MD data (circles) and from 
the exact solution of the FP equation (Eqs. [12] and [T3|) 
employing constant bulk (broken line) or true diffusivity 
profile (solid line). Agreement between solid lines and 
MD data is quite good, although deviations close to re- 
gion A are discernible and might point to residual barrier- 
crossing events orthogonal to the coordinate R Q. 
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FIG. 5. [Color online] (a) Free energy and diffusivity profiles 
in the range of separations between the first and second coor- 
dination shell for T = 300 K. (b) Commitment probabilities 
TTA and ttb- (c) Transition path probability P(TP|J?). Sim- 
ulation data (circles) are compared to FP estimates based 
on the diffusivity profile D(R) (solid lines) and a constant 
diffusivity 2_Dh20 (broken lines). 



For a one-dimensional diffusive system, the probability 
P(TP|_R) that a path passing through i? is a transition 
path (TP), i.e., a path directly connecting the regions A 
and B, is given by [l^ 

P(TP|i?) = 2^A(i?)^B(i?) = 2^A(i?)(l - 7rA(-R)), (14) 

where the factor 2 takes into account that a TP can start 
in A and reach B or vice versa. The probability P(TP|i?) 
reaches its maximum value 0.5 at the transition state de- 
noted by i?* [39^, where 7rA(i?*) = 7rB(i?*) = 0.5. A 
"good" RC is characterized by a maximum value of the 
TP probability near this one-dimensional diffusive limit 
of 0.5 [1^. In contrast, for "poor" RCs, which do not 
single out the transition states, this maximum is consid- 
erably lower; the reason is that for "poor" RCs excursions 
from and to A and excursions from and to B dominate all 
along the coordinate such that TPs are rare everywhere 
between A and B [i,[i3- 

Committor and TP probabilities are estimated by 
analyzing all simulation paths within the region R e 
[Rx^Rb] = [0.275 nm, 0.47 nm] within a time window 
of 100 ps (time resolution 5t = 0.01 ps). P(TP|i?) 
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from MD data in Fig. [5]:) reaches a maximal value of 
P{TP\Ri) « 0.38, the position Ri being slightly dis- 
placed from the FP prediction employing the D{R)- 
profile (solid line) and a constantant D (broken line) by 
about 0.02 nm; though caution is recommended in in- 
terpreting the TP probability test [13], we conclude that 
the separation R is an acceptable RC unlike in the similar 
problem of ion unbinding [9| . 



of the path perpendicular to R. Since only the magni- 
tude of this perpendicular component can be accessed, 
while the direction remains uncertain, the path R can 
not be completely reconstructed from the knowledge of 
D{R). Defining 



Ri 





dR^ 




/ dR' 


-f 




dR' 


jRa 



dR' 



D{R') 



(18) 



C. Interpretation of Diffusivity Profiles 

Based on these findings, it is possible to give a quite in- 
tuitive interpretation of the diffusivity profile. Assuming 
that the relative dynamics of two water molecules can 
be described as a diffusive process along a single path 
R in the full high-dimensional configuration space, the 
projection onto one single coordinate, in this case the 
oxygen-oxygen separation R, generally leads to consider- 
able changes in free-energy and diffusivity. 

For convenience we assume the path -R(s), s G [0, /] of 
total contour length I being arc- length parametrized, i.e., 
|di2(s)/ds| — 1 Vs, and the reactive fiux tube bing quite 
narrow such that the idea of a single, dominating path 
makes sense [ll|. The vector R = (R, R±) is split up into 
the coordinate R and an ortogonal, vectorial component 
R±, implying 



dR 




dR^ 


2 




ds 




ds 


+ 





(15) 



We assume a one-to-one correspondence between the arc- 
length variable s and the relative separation R, i.e., a 
path which does not take any value of R more than once; 
in this case, the coordinate R is just a reparametrization 
of s. As is well-known [U,!!^, such a reparametrization 
can sensibly alter free-energy and diffusivity; more pre- 
cisely the corresponding profiles along the coordinates R 
and s are connected via 



FiR) ^ Fis) + kBT log 



D{R) = D{s) 



dR 
d7 



(16) 



Combining Egs. 1151 and 1161 one deduces 



dR, 



= \ 1- 



dR 
d7 



= Wi- 



D{R{s)) 



(17) 



meaning that the knowledge of the diffusivity profile 
D{R) along the coordinate R allows to draw conclusions 
on the shape of the path R{s); that is, a reduction of the 
diffusivity D{R) along a chosen RC i? is a signature of 
pronounced contributions to the reaction path that are 
orthogonal to the RC. Deviations of the diffusivity from 
the value D{s) thus indicate a non-negligible component 



ans using the relations in Egs. [TBI and [T71 however allows 
to visualize the path in the (i?, _R^)-plane. Choosing 
a constant D(s) — 2Dii^o for illustrative purposes, we 
show in Fig.[TJ:) a fictitious path in the plane {R, R±) that 
would be consistent with the diffusivity profile D{R) ac- 
tually extracted from MD simulations. We observe that 
the pictorial reaction path has large contribution orthog- 
onal to R within the first coordination shell, where the 
diffusivity profile shows its prominent drop, i.e., for rel- 
ative separations 0.26 nm ^ R ^ 0.34 nm. Previous 
simulation results suggest that the orthogonal degree of 
freedom involved in water-bond breakage is in fact of an- 
gular nature [3, 0] ■ 



VI. CONCLUSIONS 

Summarizing, we have resolved the diffusivity profile 
for relative water dynamics based on a MFPT analysis 
of the stochastic trajectories obtained from MD simula- 
tions of SPC/E water. Over a wide range of tempera- 
tures, the diffusivity within the first coordination shell 
drops by a factor of more than six compared to large 
separations, where both partners diffuse independently. 
The form of the diffusivity profile is necessary to repro- 
duce dynamic properties such as FPT distributions ob- 
served in the simulations and can be interpreted in terms 
of a reaction path which is distorted with repect to the 
resolved separation coordinate R. 

We cautiously remark that a distorted reaction path is 
only one of a few mechanisms that would modify the lo- 
cal diffusivity; orthogonal energetic barriers (which are, 
based on our results shown in Fig. [Sj presumably small 
in the present case but dominate in related problems ^ ) 
and competin g r eaction paths (l^ or fiux-tube width 
variations (ill l4l| are additional complications. An un- 
derstanding of the precise mechanisms at work when 
water molecules separate from each other is desirable, 
but requires the applications of more complex concepts, 
which we have not pursued in this paper; among others, 
transition path theory II IJ or Markov-state modeling in 
full configuration space |42| may prove useful. The charm 
of our approach is that it allows for a consistent descrip- 
tion of the kinetics of water-bond breakage even without 
a detailed knowledge of the transition path and the in- 
volved degrees of freedom. 



8 



10-' 

fT" 10^ 

B 

10° 
10-1 

I 

+r 10"^ 

sr; 10-'' 

10"" 



1 1 1 

;_ T=280 K 


1 1 ' 


1 1 


T=300 K 






--- T=320 K 






--■ T=340 K 






'_ T— 360 K 


^^^^^^^^^^^^^^ 










" °^ ^^/{^ ^^^^ 






w 

' 1 1 


1 1 ■ 


1 ■ 



10" 



10" 



IQi 
t [ps] 



lO'' 



10-^ 



10* 



FIG. Al. [Color online] Single SPC/E water molecule MSDs 
for various temperatures {dashed color lines) and best linear 
fits to the data {solid black lines). 
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FIG. A2. [Color online] Arrhenius plot of the temperature 
dependence of the SPC/E diffusion coefficient: Symbols de- 
note results obtained through fits (see text) to the MSD data 
shown in Fig. lAf 1 the line shows that within the studied 
range of temperatures this dependence is well approximated 
by Du2o{T) ^ 956 exp (-1783.6 K/T) nm^ns. 
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Appendix A: Single Water Molecule Diffusion 
Coefficients 



Fig. ET] shows the MSD of single SPC/E molecules ex- 
tracted from MD simulations at various temperatures; 
while the MSDs show a quadratic dependence on time 
characteristic for ballistic motion on the femtosecond 
time scale, a smooth crossover to a diffusive scaling is 
observed for t « 0(ps). The diffusion constant of a sin- 
gle water molecule 



D 



H2O 



lim 



((r(t)-r(0))2) 



6t 



(Al) 



is obtained through linear fits to the data in the time 
range 10 ps < f < 10^ ps. In Tab. U our results for 
the SPC/E diffusion coefficients are compared to results 
from experiments and other simulation studies. As can 
be seen in Fig. IA2| the diffusion coefficient shows an 
Arrhenius-like temperature dependence within the inves- 
tigated range of temperatures in agreement with experi- 
mental findings [35| . 



Appendix B: Mean First-Passage Times from 
Trajectories with Finite Time Resolution 



In Fig. IB II we show several MFPT curves, which were 
obtained from the same simulation run by only varying 
the time resolution St employed for the MFPT analysis 
described in Sec. IIIIBI The reason for the differences be- 
tween the MFPT curves, are excursions beyond i?t and 
back, which are not registered due to the finite time res- 
olution St and thereby affect the estimate of the mean. 
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FIG. Bl. [Color online] Dependence of the MFPT curves 
Tfp on the time resolution St of the underlying trajectories. 
Results for several target separations Rt = 1.2 nm {orange), 
Rt = 0.8 nm {green) and Rt = 0.4 nm {blue) determined from 
MD simulation data at T = 300 K are shown. 



Note that the curves in Fig. IBll are mainly shifted ver- 
tically with respect to each other; we found that within 
the statistical uncertainty the choice of a specific time 
resolution had no visible effect on the form of the re- 
sulting diffusivity profile, which according to Eq. |S] only 
depend on the slope of the MFPT curves (comparison 
not shown). 

Note that the vertical shifts in Fig. IBll are in fact 
much larger than the resolution St of the trajectories. To 
demonstrate the influence of St on the observed MFPT 
curves, we use a simple model system: free diffusion of 
a particle along the coordinate x with diffusion constant 
Dp. A reflecting boundary at x = restricts the particle 
position to the positive part of the coordinate axis; the 
corresponding Green's function is given by 



G{x\xo;t) 



1 



4Dpt 0(3,)^ 



(Bl) 



9 



Since the process is Markovian, the probability distribu- 
tion of finding the particle at position x at time t + St is 
related to the probability distribution at time t by 



P{x;t + St) 



dx' G{x\x';St)P{x';t). 



(B2) 



According to Eq. [7]the MFPT of reaching xt when start- 
ing out from X for a flat diffusivity D{x) — Dp and free- 
energy profiles F{x) = const, is given by 



Tfp(a;,xt) = Td 1 - 




(B3) 



i/{2Dp) for diffusion 



where the characteristic time Td 
over the length xt was defined. 

We show in Fig. IB2I that MFPTs obtained from a tra- 
jectory with finite time resolution differ from the MFPT 
calculated in the continuum, because paths will eventu- 
ally cross the target and come back many times before 
a position x > Xt is first recorded in the trajectory with 
finite time resolution. A lower bound for MFPT esti- 
mates based on trajectories with finite time resolution is 
obtained by the following numeric procedure: 

1. At time < = the particle is located at xq, thus 
the probability distribution reads P{x;t = 0) 
6{x — Xq)- Consequently, the probability of finding 
the particle left of the target is Pieft — 1- Since no 
transitions beyond the target have been observed 
yet, the current MFPT estimate is ftp — 0. The 
probability distribution at time t = 5t according to 
Eg. IB2l thus simply is P{x) = G{x\xQ;5t) which is 
evaluated along x with a resolution 5x = 0.03 Xt- 

The following steps are repeated until the exit con- 
dition is met: 

2. Linearly interpolate the discrete values of P{x) to 
obtain a continous function P{x). The probabil- 
ity of still finding the particle left of the target is 
determined numerically by integration 



r>new 
Meft 



dx'P{x). 



(B4) 



3. 



4. 



The fraction x — ^icft — ^icft" particles is thus 
found on the right of the target {x > xt) for the 
first time. 

The fraction x niust have crossed the target be- 
tween time t — St and time t and contributes to 
the observed MFPT, which is updated accordingly: 
ffp = ffp -I- xit - St). 

If < 0-001, i.e., if less than 0.1 % of the par- 

ticles have not been observed right of the target at 
least once, then exit the loop and return the MFPT 
estimate ftp. Else: 



2.5 

2.0 

. ^. 1.5 

ft 
It- 

1.0 
0.5 
0.0 



1 


1 1 

St 


= 1 - 






= 0.3 


- 


St 


=0.1rd 




St 


= 0.03 Tj 




■ St 


= 


1 


1 1 





0.2 0.4 0.6 0.8 1.0 

x/xt 

FIG. B2. [Color online] MFPT curves for reaching xt in the 
case of free diffusion next to a reflecting wall at a; = for dif- 
ferent time resolution 5t of the underlying trajectory (see text 
for details). Solid colored lines are obtained by the numeric 
procedure described in the text in App. [B] the dashed black 
line denotes continuum MFPTs (Eq. IB3|I . Times are given in 
units of the characteristic diffusion time Td = a;?/(2Dp). 



(a) Numerically calculate the probability distri- 
bution at time t + St using Eq. IB2I on a dis- 
crete grid with Sx = 0.03 xt; herefore, start off 
with the interpolated version P{x; t) = P{x) 
and choose Xt as upper integration limit in 
Eq. IB2I thereby neglecting the fraction of par- 
ticles which reached separations a; > Xt in the 
last iteration. 

(b) Set Pioft = Pi"™ andt = t + St. 

(c) Go back to (2.). 



MFPT curves resulting from the procedure described 
above are shown in Fig. IB2I they show the same charac- 
teristics as the MFPT curves from MD simulation data 
in Fig. IBll namely increasing St shifts up vertically the 
curves; distortions of the curves are only observed for 
St > Tfp . As is clearly seen, one has ftp > Ttp + St for the 
smaller values of St: the deviations are thus not caused 
because the first passage time is recorded with an uncer- 
tainty on the order of the time resolution, but because 
the first observed passage time in the discretely sampled 
trajectory can exceed by far the FPT in the continous 
trajectory. 



Appendix C: Numerical Solution of the 
Fokker-Planck Equation 

When discretizing the spatial coordinate R into A'' bins 
of width AP, the FP equation (Eq. [5]) takes the form of 
a master equation (isj 



^P^{t) ^ 

dt 



W,,,_iP,_i(t) + W,,,+iP,+i(t) -f VF.,.P.(i), 



(CI) 



10 



Bin indices are denoted by z £ {1,2, . . . , N}, and the 
probability Pi{t) of observing a relative separation R 
within bin i, i.e., a separation R G [Ri — AR/2,Ri + 
AR/2), at time t as well as the transition rates Wij 
from bin j to bin i depend on both free-energy F and 
diffusivity D 



2{ARY 



2k^T 



(C2) 



where Fi = F{Ri) and 1?^ = D{Ri). Due to detailed 
balance the transition rates between neighboring bins are 
not independent 

= exp (^-^^±^^^^ VK,,,+i. (C3) 

The linear transformation Pi{t) — exp (/3Fj/2)Pi(<) with 
j3 = {k-QT)^^ converts Eg. ICll into a simple differential 
equation involving a tridiagonal, symmetric matrix with 
entries W,, 



^P^{t) 
dt 



N 



^ W^j Pj{t) with W,j = e'^^' /2 W^J 



PF,/2 



(C4) 

which is readily solved in terms of a matrix exponential 

N 

P,(<) = ye-^^'/2(e^*) e^^^/2p,(0) (C5) 

JV 

^^G.,(i)P,-(0), 

i=i 

where in the last equality the Green's function Gij speci- 
fying the probability of landing in bin i a time t after the 
given start in bin j was defined. The matrix exponential 
in Eq. IC5I is computed numerically by diagonalization of 
the symmetric matrix W = QAQ^^, with Q being the 
matrix of eigenvectors and A being the diagonal matrix 
of eigenvalues of W. The matrix exponential is then sim- 
ply given by 



Qe^*Q-\ (e^*) =5,,e 



(C6) 



For the case of relative SPC/E water dynamics, a bin 
width AR — 0.002 nm was used; a reflective boundary 
condition (Wo,i — W^i.o — 0) was imposed at i?niin = 
0.235 nm corresponding to a value of the free-energy F « 
18 kBT. 

FPT distributions are obtained by imposing an absorb- 
ing boundary condition at the target position Rt, thus 
disregarding paths in the time evolution of Gij which 
have already reached the target beforehand: the total 
number N of bins is chosen such that the target Rt is 
part of bin N + 1 and the absorbing boundary condition 
is implemented by setting 



Wn,n — —Wn-i,n — W^Af+l, 



N, 



(C7) 



but neglecting the backward flux WN_N+iPN+i{t) since 
PN+i{t) = Vt. The survival probability Pi^,.^ for a start 
in bin j is 



N N 
1=1 i=l ' 



(C8) 

which is evaluated numerically at times t G [0, 200 ps] 
with time resolution St = 0.1 ps. The first-passage time 
(FPT) distribution is approximated by the finite differ- 
ence 

r /, I r , /o ■^ ^ ^s'urv(^) " Piuivi^ ^ fr<r\\ 

}ip[t + bt/2;]) ^ — . (C9) 

The MFPT is obtained by taking the first moment of the 
FPT distribution, 



poo poo 

Tfp= / dttffp{t)= / dtPsurvW 
Jo JQ 



M-1 



St 



(CIO) 



where in our case M = 2000. 

FPT histograms from MD data and from the numerical 
solution of the FP equation described in this section for 
other target separations than the one shown in Fig. HJd) 
are found in Fig. ICll No significant impact on the FPT 
distributions was observed when refining the discretiza- 
tion in space and/or time. 
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